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1  Introduction 


The  main  objective  of  this  project  is  the  development  of  accurate  physical 
models  and  efficient  numerical  algorithms  suitable  for  diffraction  tomographic 
reconstruction  of  the  compressibility,  acoustic  attenuation,  and  mass  density 
of  the  prostate  through  a  multielement  transurethral  ultrasound  transceiver. 
Unlike  conventional  ultrasound  imaging,  which  is  non-quantitative  and  af¬ 
fected  by  speckle  artifacts,  three-dimensional  tomographic  solution  of  the 
inverse  acoustic  scattering  problem  has  the  potential  to  quantitatively  recon¬ 
struct  the  detailed  acoustic  properties  of  the  prostate  from  measurements  of 
the  scattered  radiation  field.  However,  numerous  challenges  must  be  con¬ 
fronted  before  such  a  proposition  becomes  feasible  in  real-world  applications. 
It  is  well-known  that  problems  in  inverse  scattering,  in  addition  to  being 
mathematically  complex  and  computationally  intensive,  are  also  particularly 
ill-posed  and  ill-conditioned.  [1,  2,  3]  This  progress  report  takes  up  where  our 
previous  report,  for  the  award  period  from  March  1,  2004  through  Febru¬ 
ary  28,  2005,  left  off,  so  we  minimize  duplication  of  results  already  reported 
previously. 

Our  research  over  the  past  year  has  focused  on  three  primary  goals:  1) 
implementation  and  testing  of  various  schemes  for  solving  inverse  scattering 
problems  in  2D  and  3D  and  their  optimization;  2)  development  and  test¬ 
ing  of  a  theoretical  foundation  and  practical  implementation  of  the  acoustic 
scattering  problem  in  the  endoluminal  geometry  in  which  the  ultrasound 
transceivers  are  located  within  the  lumen  of  the  urethra  (which  we  term  en¬ 
doluminal  ultrasound  tomography ,  ELUST),  and  3)  development  of  analytical 
and  numerical  contrast-source  inversion  procedures  for  simultaneous  recovery 
of  sound  speed  and  mass  density.  [4,  5,  6,  7,  8,  9] 

2  Inversion  in  2D  and  3D 

2.1  Inverse  Modeling 

We  have  developed  and  implemented  a  nonlinear  conjugate  gradient  (NLCG) 
algorithm  [10,  11]  for  minimization  of  a  two-term  Tikhonov  cost  functional  [3] 
incorporating  both  a  y2  term  for  model-data  agreement  and  a  regularization 
term  for  object  error.  The  advantage  of  NLCG  algorithms  over  quasi-Newton 
methods  is  that  they  only  require  the  evaluation  of  the  gradient  of  the  ob- 
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jective  function,  rather  than  its  Hessian,  potentially  leading  to  significant 
computational  and  storage  savings  [10]  in  cases  such  as  ours  where  function 
evaluations  are  expensive.  We  further  improve  the  efficiency  of  our  minimiza¬ 
tion  scheme  by  implementing  a  gradient  calculation  based  on  the  method  of 
adjoint  fields  that  allows  us  to  forego  the  calculation  of  the  full  Jacobian 
matrix  at  each  iteration.  In  our  work  to  date  we  have  considered  both  L-2- 
norm  and  total  variation  based  Li-functionals.  The  former  tends  to  enforce 
a  higher  degree  of  smoothness  on  the  solutions,  while  the  latter  is  known  to 
have  better  edge-preserving  properties  and  is  well  suited  to  reconstruction 
of  block  objects  with  sharp  boundaries.  Both  methods  provide  reasonably 
good  results  in  the  test  phantoms  we  have  considered,  although  total  vari¬ 
ation  appears  to  be  somewhat  superior  for  the  limited  test  cases  studied  so 
far.  More  work  remains  to  be  done  to  establish  definitively  if  one  is  superior 
to  the  other  for  ELUST  applications. 

Care  must  be  taken  in  the  testing  of  the  inversion  algorithms  to  not  com¬ 
mit  an  “inverse  crime”. [1,  2]  This  happens  when  the  same  forward  model  is 
used  to  compute  both  “measured”  and  reconstructed  scattered  fields.  In  this 
case,  it  is  possible  to  obtain  spuriously  good  agreement  of  the  reconstructed 
and  true  objects  and  be  misled  into  overestimating  the  accuracy  of  the  in¬ 
version  scheme.  For  single  or  multifrequency  reconstructions  at  a  fixed  grid 
size,  we  perform  the  following  steps  in  computing  the  reconstructed  object: 

•  Calculate  simulated  measurement  data  from  “true”  object. 

The  “true”  object  is  discretized  on  a  finer  grid  than  the  highest  resolu¬ 
tion  reconstruction  desired  (in  our  simulations,  typically  twice  Nrecon ). 
For  example,  if  we  wish  to  reconstruct  a  phantom  on  a  64x64  grid,  the 
measured  scattered  held  is  computed  from  an  discretized  phantom  that 
is  no  smaller  than  128x128. 

•  Add  noise  to  measurement  data:  Gaussian-distributed  random 
noise  with  a  standard  deviation  of  ~  2  —  3%  is  separately  added  to  the 
real  and  imaginary  components  of  the  measured  scattered  held. 

•  Determine  reconstructed  object  from  inversion  algorithm:  The 

cost  functional  is  minimized  using  the  NLCG  routine  to  give  the  recon¬ 
structed  object. 
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2.2  Multigrid-Multifrequency  Inversion 

A  major  impediment  to  the  development  of  effective  inverse  scattering  algo¬ 
rithms  is  the  tendency  of  the  solutions  to  become  trapped  in  local  minima 
of  the  objective  function,  leading  to  reconstructions  that  differ  significantly 
from  the  true  object  and  dependence  on  the  starting  guess.  In  order  to  com¬ 
bat  this  problem,  we  have  developed  a  multigrid,  multifrequency  technique 
that  appears  to  work  quite  effectively.  In  some  cases,  the  effective  regulariza¬ 
tion  implicit  in  iteratively  bootstrapping  from  low  to  high-resolution  appears 
to  obviate  the  need  for  a  regularization  term  entirely.  This  technique  is 
implemented  as  follows: 

•  Choose  the  starting  grid  size  for  an  initial,  low-resolution  re¬ 
construction.  Typically,  N0  ~  16. 

•  Choose  a  uniform  object  for  the  starting  guess,  70.  Typically, 
70  —  0.01. 

•  k-th  multigrid  iteration: 

—  Set  wavelength  for  reconstruction:  Typically,  A k/L  =  4 /Nk, 
where  L  is  the  physical  dimension  of  the  object,  and  A&,  Nk  are 
the  wavelength  and  grid  size,  respectively,  in  the  k-th  iteration. 
Such  a  choice  of  A*,  and  Nk  allows  the  Nyquist  sampling  interval 
to  be  exceeded  and  thereby  ensuring  that  the  forward  model  is 
able  to  accurately  represent  the  scattered  held. 

—  Perform  reconstruction:  7*,  is  reconstructed  by  FFT-NLCG. 

—  Increment  grid  size:  Nk+1  =  (3Nk.  Typically,  j3  —  1.25  —  2.00. 

—  Compute  an  improved  starting  guess:  Upsample  7*.  to  reso¬ 
lution  Afc+i  to  obtain  jk+ 1  • 

•  Terminate  iteration:  Finished  when  Ay.  reaches  the  ultimate  desired 
size  for  the  reconstruction. 

In  multigrid  reconstructions,  the  measurement  data  for  each  resolution  (and 
the  corresponding  frequency)  is  computed  on  the  same  grid  size,  normally 
twice  the  largest  resolution  to  be  reconstructed.  This  ensures  that  the  inver¬ 
sion  is  not  affected  by  discretization  artifacts  that  change  with  lattice  size. 
The  effectiveness  of  the  multigrid  method  is  demonstrated  in  Fig.  1,  where 
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we  show  the  reconstructions  of  a  Shepp-Logan  phantom  for  a  grid  size  of 
128x128  using  this  approach.  In  this  case,  grid  dimensions  of  16,  20,  24,  28, 
32,  40,  48,  56,  64,  and  128  were  used  sequentially,  and  each  reconstruction 
was  computed  using  three  wavelengths  simultaneously:  A  =  L,  2 L\fiV,  and 
4 L/N.  In  the  figure,  panel  (a)  shows  the  true  object  at  256x256  grid  res¬ 
olution,  and  panel  (b)  shows  the  reconstructed  object.  It  is  apparent  that 
we  are  able  to  resolve  objects  down  to  a  size  of  approximately  4  pixels,  as 
expected  from  the  minimum  wavelength  used. 

Multigrid  inversion  is  compared  with  direct  inversion  without  multireso¬ 
lution  bootstrapping  in  Figure  2.  The  upper  panels  (a-c)  are  reconstructions 
of  the  scatterer  shown  in  Figure  3d  performed  without  regularization  (a  =  0) 
with  a  uniform  starting  guess  for  three  different  grid  sizes  and  wavelengths 
defined  by  A /L  —  4L/N.  Here,  the  Gibbs  phenomenon  is  highly  apparent, 
with  the  actual  object  ringed  with  large  amplitude  haloes.  The  lower  panels 
(d-f)  are  identical  except  that  multigrid  iteration  was  utilized  so  that  the 
starting  guess  for  Figure  2e  was  an  upsampled  version  of  the  reconstruction 
in  Fig.  2d  and  the  starting  guess  for  2f  was  an  upsampled  version  of  3e. 
Here  the  regularizing  effect  of  multigrid  iteration  is  clearly  apparent.  The 
reconstructions  in  Figure  2a  and  Figure  2d  are  identical,  but  use  of  the  lower 
frequency  reconstructions  as  an  initializer  for  the  higher  frequencies  leads  to 
dramatic  improvement  in  the  agreement  between  the  true  scatterer  and  the 
reconstructed  objects. 

2.3  Inversion  in  the  Exterior  Source/Interior  Detector 
Geometry 

In  this  geometry,  the  source  is  a  plane  wave,  as  previously,  but  the  detectors 
are  embedded  in  the  object  to  be  reconstructed.  If  it  is  assumed  that  the 
detectors  are  pointlike  objects,  and  do  not  significantly  perturb  the  scattered 
fields,  then  the  forward  problem  can  be  treated  in  a  manner  similar  to  the 
conventional  exterior  source/exterior  detector  (or  exterior-exterior)  geome¬ 
try.  The  main  difference  is  that  the  scattered  fields  are  now  the  same  as  the 
total  internal  fields. 

In  order  to  perform  a  feasibility  test  for  reconstructing  in  an  endoluminal 
geometry,  we  implemented  a  mixed  geometry  simulation  in  which  plane  wave 
sources  are  incident  on  the  object  from  outside,  but  the  detectors  (which 
are  regarded  as  pointlike  and  do  not  perturb  the  scatterer)  are  arrayed  in 
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the  interior  of  the  object.  As  pointed  out  above,  in  the  exterior /interior 
geometry,  the  internal  fields  are  to  be  used  for  inversions  rather  than  the 
scattered  fields.  We  have  considered  both  a  standard  modified  Shepp-Logan 
phantom  (Fig.  3a-c),  and  a  phantom  intended  to  more  closely  resemble  a 
real  endoluminal  configuration  (panels  d-f).  In  both  cases,  we  note  that, 
although  the  effective  maximum  frequencies  are  relatively  low  (~2  MHz), 
the  reconstructions  clearly  show  the  overall  structure  of  the  object  and,  in 
addition,  reproduce  many  of  the  fine  structures  quite  well. 

2.4  Anisotropic  Diffusion  Filtering 

As  a  consequence  of  the  finite  number  of  wavelengths  used  in  reconstruction, 
images  typically  manifest  a  type  of  Gibbs  ringing.  Using  multifrequency  data 
for  inversion,  as  in  the  reconstructions  in  Fig.  1,  mitigates  this  significantly, 
though  at  the  price  of  an  increase  in  the  computation  time.  The  ringing  is 
particularly  apparent  in  panel  (b)  of  the  single  frequency  reconstruction  of 
Fig.  3,  where  it  is  visible  both  within  and  outside  the  reconstructed  object. 
A  method  for  reducing  noise  in  multidimensional  data  while  simultaneously 
preserving  sharp  edges  which  has  come  into  relatively  wide  use  in  image 
processing  and  computer  vision  in  the  past  decade  is  known  as  anisotropic 
diffusion  filtering  (ADF)  [12,  13].  Because  we  have  a  reasonable  physical 
estimate  of  the  length  scale  of  these  spurious  oscillations,  ADF  is  particularly 
easy  to  apply  as  a  post-processing  step  once  reconstruction  is  complete.  The 
ability  of  ADF  to  remove  the  high-frequency  noise  without  compromising  the 
representation  of  the  thin,  strongly  scattering  shell  is  clearly  demonstrated  in 
panel  (c)  of  Fig.  1  as  well  as  panels  (c)  and  (f)  of  Fig.  3.  We  have  investigated 
the  iteration  properties  of  ADF  of  various  reconstructions  and  find  that,  as 
shown  in  Figure  4,  there  is  typically  a  minimum  in  the  RMS  error  between  the 
true  and  reconstructed  objects  for  a  finite  number  of  ADF  iterations  in  the 
range  of  5-15.  It  appears  that  integration  of  a  cycle  of  ADF  into  the  multigrid 
reconstruction  algorithm  minimizes  the  appearance  of  spurious  oscillations 
and  leads  to  enhanced  reconstruction  accuracy,  though  future  work  is  needed 
to  establish  the  optimal  parameter  choices. 

2.5  Three-dimensional  Reconstruction 

Figures  5-7  demonstrate  the  capability  to  perform  full  3D  reconstructions  us¬ 
ing  the  Lippmann-Schwinger  forward  model.  In  this  case,  the  reconstruction 
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was  performed  in  the  exterior-exterior  geometry,  with  256  detectors  arrayed 
in  4  rings  of  64  lying  perpendicular  to  the  z-axis.  While  the  computational 
demands  limited  the  size  of  our  reconstruction  to  24  x  24  x  24,  so  discretiza¬ 
tion  and  the  relatively  long  wavelengths  used  make  conditions  sub-optimal 
for  inversion,  Figure  ??  shows  the  ability  to  accurately  depict  the  large  scale 
interior  structures  of  the  object  accurately,  and  appears  to  perform  as  well 
as  the  2D  algorithm  at  modeling  the  thin  exterior  shell.  Other  than  reformu¬ 
lation  of  the  scattering  Green’s  function  and  use  of  three  dimensional  source 
and  detector  arrays,  the  3D  problem  is  handled  by  the  same  inversion  code 
as  2D  problems. 


3  Theoretical  Developments 

3.1  Scattering  in  the  Interior-Interior  Geometry 

The  interior- interior  scattering  geometry  is  a  novel  geometry  that  is  not  con¬ 
ventionally  considered  in  inverse  scattering  problems.  In  most  cases,  the 
source  as  well  as  detectors  are  exterior  to  the  scatterer.  In  the  endoluminal 
tomography  (ELUST)  of  the  prostate,  sources  and  detectors  are  both  internal 
to  the  organ.  While  our  preliminary  work  will  focus  on  extension  of  existing 
algorithms  to  the  case  of  nonperturbing  interior  point  sources,  application  of 
this  methodology  in  an  experimental  setting  will  require  treatment  of  the  per¬ 
turbation  of  the  acoustic  scattering  caused  by  the  presence  of  the  catheter 
tip  and  transceiver  array.  We  have  developed  a  theoretical  foundation  for 
the  case  where  it  is  assumed  that  the  transducers  are  placed  on  the  surface 
of  a  sound-impermeable  cylinder  located  within  the  urethra.  We  treat  this 
cylinder  assuming  Neumann  boundary  conditions  where  the  normal  compo¬ 
nent  of  wave  velocity  on  the  surface  vanishes.  In  the  usual,  exterior-exterior 
geometry,  the  held  propagator  is  the  outgoing  free-space  Green’s  function, 
but  the  presence  of  the  Neumann  cylinder  changes  the  situation  in  a  compli¬ 
cated  manner.  The  occurrence  of  multiple  scattering  not  only  in  the  tissue 
matrix  alone,  but  between  the  tissue  and  the  cylinder  makes  it  necessary 
to  modify  the  free-space  Green’s  function.  The  resulting  Neumann  Green’s 
function  must  obey  the  Helmholtz  equation  and  satisfy  both  the  Neumann 
boundary  condition  and  Sommerfeld’s  radiation  condition  at  infinity.  For  an 
arbitrary  scatterer,  such  a  Green’s  function  is  difficult  to  construct,  but  in 
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the  canonical  cylindrical  geometry  it  is  possible  to  obtain  G n- 
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in  3-D.  In  Eq.  (1),  x  =  (r,  (ft),  and  y  =  (r,(j)').  Moreover,  r<(>)  implies  the 
larger  or  smaller  of  r  and  r' .  h0  is  the  incident  wavenumber,  and  J,  H  are 
the  standard  cylindrical  Bessel  and  Hankel  function,  respectively  [14].  In  Eq. 
(2),  x  =  (r,  (j),z),  y  =  (V,  0',  z1),  k±  =  k2x  +  ky ,  and  k\\  =  s/k^  —  k\.  a  is  the 
radius  of  the  cylinder.  Cylindrical  coordinates  were  used.  The  expression 
in  Eq.  (2)  is  for  the  supersonic  waves  only.  For  subsonic  waves,  where  the 
waves  are  evanescent ,  the  cylindrical  Hankel  functions  must  be  replaced  by 
the  modified  Bessel  functions  Kn  [15]. 

The  expression  for  the  incident  held  resulting  from  an  infinitesimal  rect¬ 
angular  transducer  on  the  cylinder  surface  has  also  been  derived.  We  present 
only  the  3-D  result.  Let  the  transducer  be  2 L  long,  located  at  a  distance  of  Z 
from  the  origin  of  the  cylindrical  coordinate  system  on  the  cylinder,  and  have 
an  angular  width  of  20O--  Moreover,  let  Vo,  do,  Co,  Vo  represent  the  amplitude 
of  vibration  of  the  transducer  surface,  the  mass  density,  sound  speed  in  the 
homogeneous  ambient  medium,  respectively.  Then  the  incident,  i/>(x),  at  a 
point  x  =  (r,  (j),  z),  was  derived  as: 


x )  =  1 0±^°L  (id0c0k0)  ^  ein<t> sinc{;n(j) 0)  • 
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dk\\  sinc{k\\L ) 


n=— oo 

kb 


Gn  ( k±p )  Jkp(z-(Z+L)) 


Ir 1  Gn\k±a) 


(3) 


sinc(x )  =  sinx/x,  is  the  sinc-function.  Now,  for  all  practical  purposes,  both 
L  and  <f>o  are  expected  to  be  small.  In  that  case,  the  sinc-functions  can  be 
approximateded  by  unity,  and  Eq.  (3)  reduces  to 

ip{in\x)  =  2iv0±<j)0d0c0  (4) 

Gn  (fc0  -L  a) 

Note  that  the  held  given  in  Eq.  (4)  is  roughly  in  the  plane  passing  through  the 
midplane  of  the  transducer.  This  is  interesting  since  it  raises  the  possibility 
of  doing  3-D  imaging  via  reconstructing  only  2-D  planes,  a  much  simpler 
proposition. 

With  Gn  given  in  Eq.  (2),  and  in  Eq.  (3)  or  Eq.  (4),  the  scattered 
held  anywhere  in  the  prostate  can  be  obtained  by  the  Lippmann- Schwinger 
integral  equation  of  scattering,  which  is: 

ip(x)=ip(m)(x)  +  [  GN{x\y)'y{y)ip(y)dy.  (5) 

Jn 

4  Simultaneous  Reconstruction  of  Compress¬ 
ibility  and  Mass  Density 

In  ultrasonic  medical  diagnostics,  essentially  all  attempts  at  quantitative 
acoustic  imaging  have  focused  on  speed  of  sound  and  attenuation  imaging,  [17, 
18,  19,  20]  while  neglecting  mass  density  (p)  variation.  Addressing  the  latter 
is  relatively  more  difficult  than  handling  the  compressibility  (k),  because  of 
the  involved  manner  in  which  the  density  term  appears  in  the  basic  integral 
equation  of  scattering.  Unlike  the  compressibility,  the  density  term  involves 
differential  operators,  and  consequently,  simultaneous  determination  of  com¬ 
pressibility  and  mass  density  from  the  scattering  data  is  complicated.  [16] 
This  is  the  reason  that  the  mass  density  reconstruction,  when  undertaken, 
is  conventionally  limited  to  the  linearized  Born-Rytov  approximation  or  for 
uniform  distributions.  [21,  22,  23]  We  have  developed  a  theoretical-numerical 
procedure  for  inverting  both  compressibility  and  density  distributions  under 
non-Born  conditions. 

The  full  scattering  problem  involving  both  (p)  and  (k)  is  given  by  [24]: 

(A +  k20)p=  -kfocp  +  V  •  (7 PVp)  (6) 
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plus  Sommerfeld’s  radiation  condition,  namely 

limit\x\->00\x\m  [d\x\ps  -  ik0psi\  =  0, 

The  corresponding  integral  equation  of  scattering  is: 

Ps=  [  G(0\x\y )  [—k^cp  +  V  •  (7PVp)]  ,  x  £  ft.  (7) 

The  free-space  Green’s  function,  is  {i / A)H^\kQ\x— y\)  in  2D  and  [l/47r|a;— 
y^elk°\x-y\  in  3D.  7C  =  1  —  and  =  1  —  po/ p,  are  the  deviations  of  k 

and  p  relative  to  the  background. 

We  have  developed  the  equivalent  source  based  recent  method  of  contrast- 
source  inversion  (CSI).  An  earlier  version  of  a  CSI-likc  method  was  consid¬ 
ered  by  Colton  [1]  in  the  context  of  obstacle  scattering.  The  contrast-source 
inversion  avoids  solving  the  CPU-intensive  full  forward  problem.  Moreover, 
it  employs  only  linear  CG,  in  which  the  computationally  most  involved  step 
length  calculations  can  be  performed  analytically.  In  CSI,  the  objective  func¬ 
tional  is  constructed  as: 

$(7c7p)  =  lYl 
3= 1 

The  summation  index  j  runs  over  the  incident  waves  which  are  assumed  to  be 
J  in  number.  &Dj,  and  &oj  represent  the  j-th  least  squares  constraints  on  the 
data  and  the  object  error,  respectively.  Specifically,  =  || <&D  —  G^*(j)j\\2D, 
and  the  object  residual  is  —  ||-7r[p}n  —  G^  *  4>j  —  (j>j  is  the 

j-th  equivalent  source.  C  and  D  stand  for  the  computational  and  the  de¬ 
tector  domain,  respectively.  The  normalization  constants  in  Eq.  (8)  are: 
$£)  =  ]W7=1  I^dI2  and  $o  =  J2'j=i  \Tpinf  .  T  is  an  operator  given  by: 
T  =  [— &o7c  +  V  (7PV)] .  A  comparison  of  T  so  defined  with  its  expression: 
T  =  —  /cq7c  when  only  n  is  to  be  reconstructed,  demonstrates  the  complexity 
involved  in  inverting  mass  density.  The  iterative  inversion  proceeds  by  gener¬ 
ating  the  sequences,  {4>jn},  {7 c}n,  and  {7 p}n.  Assume  that  {4>j,n- 1},  {7c}n-i, 
and  {7P}n_i  are  known.  The  source,  <j)j  will  be  updated  first  after  which  the 
object  parameters  will  be  updated. 

We  have  analytically  calculated  the  key  parameters  needed  for  The  CSI 
minimization  of  the  functional  in  Eq.  (8)  requires  analytical  expressions  for 
some  key  parameters.  We  have 


$Dj  +  $Oj 


$ 


D 


$ 


O 


(8) 
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Developed  the  full  analytic  Frechet  derivative  of  <f>  in  (f)Jyn- 1  • 


•  Partially  completed  the  analytic  step  lengths,  a%,  a£,  and  apn  needed 
in  updating  the  contrast  and  the  source. 

•  We  have  also  completed  the  modifications  that  are  needed  in  (1)  and 
(2)  for  complex  compressibility,  i.a.,  including  attenuation. 

The  basic  calculations  for  CSI  inversion  of  (p)  and  (k)  will  be  complete 
shortly,  and  numerical  implementation  will  be  ready  to  be  undertaken. 

5  Key  Research  Accomplishments 

The  major  accomplishments  over  the  time  period  reported  are  inversions  in 
3-D  and  scattering  and  inverse  scattering  in  the  novel  semi  and  fully  endolu- 
minal  prostate  scattering  geometry.  The  key  research  accomplishments  are 
itemized  below. 

•  Numerical  inversions  in  three  dimensions. 

—  Numerical  implementation  of  3-D  Green’s  function  [25,  26]. 

—  Numerical  implementation  of  CG-FFT  [27,  28,  29]  3-D  Lippmann- 
Schwinger  forward  solution  using  Richmond’s  scheme  of  discretiza¬ 
tion  [30]. 

—  Verifications  of  3-D  forward  solutions  against  the  analytical  exact 
3-D  results. 

—  Numerical  NLCG  3-D  inversion 

•  Initial  investigation  of  various  forms  of  regularizing  functional  and  var¬ 
ious  approaches  to  optimizing  the  convergence  of  minimization  of  the 
objective  function. 

•  Minimization  of  oscillatory  reconstruction  artifacts  using  anisotropic 
diffusion  filtering. 

•  Numerical  2-D  inversions  in  the  novel  exterior /interior  semi-endoluminal 
prostate  scattering  geometry. 
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•  Theoretical  formulation  of  forward  scattering  in  the  novel  fully  endo- 
luminal  prostate  scattering  geometry. 

•  Development  of  theoretical  and  numerical  procedures  for  simultaneous 
reconstructions  of  compressibility  and  density  by  the  equivalent  source 
based  contrast-source  inversion  technique. 


6  Reportable  Outcomes 

•  Abstract  and  oral  presentation  at  the  17th  Annual  UCAIR  Symposium, 
Park  City  UT,  October  2005. 

•  Abstract  presented  at  2005  Fully  Three-Dimensional  Image  Recon¬ 
struction  Meeting  in  Radiology  and  Nuclear  Medicine,  Salt  Lake  City 
LIT,  July  2005. 

•  Invited  talk  on  the  application  of  the  method  of  adjoint  fields  to  inverse 
scattering  at  2006  Acoustical  Society  of  America  Annual  Convention, 
Providence  RI,  June  2006. 

•  Paper  on  the  analysis  of  error  propagation  in  nonlinear  diffraction  to¬ 
mography  to  be  submitted. 

7  Conclusions  and  Future  Work 

From  the  discussions  above,  it  is  apparent  that  a  sufficiently  accurate  rep¬ 
resentation  of  the  prostate  will  require  a  resolution  for  which  the  compu¬ 
tational  burden  will  exceed  reasonable  limits  for  a  single  processor  system. 
Fortunately,  much  of  the  CGFFT  computation  is  readily  parallelized  to  take 
advantage  of  multiprocessor  hardware  and/or  high-performance  computing 
cluster  systems  that  are  rapidly  growing  in  popularity  and  availability.  Nev¬ 
ertheless,  while  such  approaches  are  likely  to  be  relevant  to  this  work,  it  still 
remains  critical  to  investigate  all  avenues  to  improving  the  intrinsic  algorithm 
performance  before  resorting  to  brute  force.  For  the  upcoming  year,  we  will 
work  in  a  number  of  areas: 
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7.1  Improving  Forward  Model  Performance 

•  develop  preconditioning  strategies  for  acceleration  of  conjugate  gradient 
convergence. 

•  implement  the  non-Uniform  FFT  [27,  ?]  to  take  advantage  of  unevenly 
spaced  grids. 

•  predictor-corrector  methods  for  minimizing  evaluations  of  full  forward 
model 

•  parallelization  of  forward  model  algorithm  and  testing  on  cluster  system 

•  Further  develop  the  theory  and  implement  ELUST  geometry  for  the 
prostate. 

•  Further  develop  theory,  numerics  and  implement  simultaneous  inver¬ 
sion  of  complex  compressibility  and  mass  density. 

7.2  The  Inverse  Problem 

We  have  demonstrated  the  inversion  efficacy  of  NLCG  with  the  adjoint  held 
gradient  computation  [31]  by  implementing  it  in  the  reconstructions  of  ob¬ 
jects  in  both  two  and  three  dimensions.  However,  there  is  much  improvement 
to  be  made  in  the  inversion  algorithm.  Among  the  strategies  which  we  intend 
to  employ  over  the  next  year  are: 

•  investigation  of  alternative  and/or  hybrid  minimization  strategies  such 
as  simulated  annealing,  genetic  algorithms  and  CSI  type  techniques  to 
be  used  in  conjunction  with  nonlinear  conjugate  gradient  methods  to 
maximize  convergence  rate. 

•  investigate  the  effects  of  different  regularizing  functionals  including  to¬ 
tal  variation,  maximum  entropy,  laplacian,  and  others,  on  the  inverse 
problem  convergence  properties  of  inversion  algorithms. 

•  parallelization  of  inversion  algorithms. 

•  development  of  a  realistic  prostate  phantom  using  MRI  and  conven¬ 
tional  ultrasound  data  in  conjunction  with  anatomic  information. 
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•  implementation  and  testing  of  2.5D  reconstruction  using  3D  forward 
modeling  to  compute  the  scattered  field  in  conjunction  with  sequential 
planar  2D  reconstructions  of  the  scatterer. 
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8  Figures 
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Figure  1:  Reconstruction  in  the  exterior  geometry  of  a  2D  scattering  phantom 
on  a  128x128  lattice.  Panel  (a)  shows  the  true  object,  based  on  the  modified 
Shepp-Logan  phantom.  Panel  (b)  shows  the  object  reconstructed  by  our  nonlinear 
conjugate  gradient  inversion  algorithm  with  L 2  norm  regularization,  24  incident 
plane  waves  with  three  wavelengths,  and  an  array  of  256  detectors  located  on  a 
ring  at  a  distance  of  Rd  =  4 L  from  the  center.  Panel  (c)  shows  the  reconstruction 
from  (b)  with  8  iterations  of  anisotropic  diffusion  filtering  applied  to  reduce  Gibbs 
artifacts.  The  maximum  scattering  inhomogeneity  of  the  object  is  7  =  0.3  and  all 
three  images  are  plotted  on  an  identical  intensity  scale. 
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Figure  2:  Comparison  of  reconstructions  in  the  exterior/interior  geometry  with¬ 
out  (panels  a-c)  and  with  (panels  d-f)  multigrid  iteration,  clearly  demonstrating 
the  regularizing  effect  of  this  approach.  The  true  object,  shown  in  panel  (d)  of 
Fig.  3,  was  discretized  on  a  128x128  grid.  Reconstruction  was  performed  with  64 
incident  waves  and  64  detectors  located  as  in  the  previous  figure;  the  regularization 
parameter,  a,  was  set  to  zero  in  both  cases. 
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Figure  3:  Reconstructions  in  the  exterior/interior  geometry.  Panels  (a-c)  show  the 
56x56  phantom  with  7 mai  =  0.3  (panel  a),  the  multigrid  reconstruction  for  128 
incident  plane  waves,  and  an  array  of  64  detectors  located  on  a  ring  at  a  distance 
of  Re>  =  L/8  from  the  center  of  the  object  (panel  b),  and  the  reconstruction 
from  (panel  b)  after  application  of  5  iterations  of  anisotropic  diffusion  filtering  to 
reduce  Gibbs  artifacts.  Panels  (d-e)  show  the  same  for  a  48x48  phantom  more 
representative  of  the  endolunrinal  geometry.  The  reconstruction  in  panel  (b)  was 
performed  with  64  incident  plane  waves  and  64  detectors  in  the  same  configuration 
as  panels  (a-c).  All  objects  are  plotted  on  an  identical  intensity  scale. 
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RMS  Error 


1.75 


Figure  4:  Residual  root  mean  square  (RMS)  error  between  reconstructed 
and  true  object  for  the  reconstruction  shown  in  Figure  1  as  a  function  of 
iterations  of  anisotropic  diffusion  filtering.  The  presence  of  a  clear  minimum 
for  a  finite  number  of  iterations  (in  this  case  around  11)  indicates  that 
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Figure  5:  Sections  through  a  24  x  24  x  24  numerical  phantom  developed 
as  a  3D  analog  to  the  2D  acoustic  Shepp-Logan  phantom  for  testing  three- 
dimensional  reconstruction  algorithms.  Here,  the  maximum  value  of  7  was 
0.3. 
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Figure  6:  Sections  of  a  multigrid  reconstruction  in  3D  of  the  phantom  shown 
in  Figure  5.  Here  the  total  variation  regularization  term  was  used  with 
a  =  1CT3.  The  scale  is  identical  to  that  of  Figure  5. 
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Figure  7:  Reconstruction  error  for  the  3D  phantom  plotted  in  Figure  5, 
plotted  from  -0.06  to  0.06. 
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